function xdot = rocketODE(t,x,data)
% function to return derivatives for water bottle rocket state equations
% inputs are t = time (scalar), x = state vector (see below), and data
% (rocket adn atmospheric conditions)
%
% necessary components in data structure are:
% data.tt, data.thrust -> time and thrust vectors from static test data,
% interpolated, if no thrust data (solely glide, such as necessary for ISP
% approximations, initialize with tt = [0,1000], thrust = [0,0]
%
% data.wind -> 3-component vector of wind velocity
% data.initH -> normalized vector of initial heading when within launch pad
% data.bottleMass -> mass of empty rocket
% data.airDens -> density of air (atmosphere)
% data.bottleXArea -> rocket cross-sectional area for drag computation
% data.CD -> rocket Coefficient of drag
% data.gavec -> 3-component vector of gravitational acceleration 
% data.rhoW -> density of water or 'rocket fuel'....
% data.nozzleArea -> cross sectional area of the nozzle/bottle neck/exit


% x vec is [x,y,z,vx,vy,vz,mass water]
 
% magnitude of position vector
posMag = norm(x(1:3));

%calculate thrust at current time
thrust = linterp(data.tt,data.thrust,t);

% calc booleen as to whether rocket is still within pad (within 1m)
withinPad = posMag < 1;

% calculate relative velocity (take into account wind)
vrel = x(4:6) - data.wind;

%magnitude of wind relative velocity
vrelMag = norm(vrel);

%calculate heading, uses initial heading if within pad and the relative
%current velocity if not within pad
h = data.initH*withinPad + ~withinPad*vrel/(vrelMag + 1e-13);

totalMass = data.bottleMass + x(7); %total mass of rocket is the bottle mass in addition to the water mass

%calculate drag magnitude
Drag = 0.5*data.airDens*(vrelMag^2)*data.bottleXArea*data.CD;

%calculate gravity force
gvec = data.gavec*totalMass;

%caculate water mass flow
vwater = sqrt(thrust/(data.coeffDischarge*data.rhoW*data.nozzleArea));
mdotw = data.coeffDischarge*data.rhoW*data.nozzleArea*vwater;

%check to see if there's any water left
isEmpty = x(7) <= 0;

%set mass flow of water equal to zero if bottle is empty
mdot = mdotw*~isEmpty;


%% lets output stuff
xdot(1:3) = x(4:6);

xdot(4:6) = (thrust*h - Drag*h + gvec)/totalMass;

xdot(7) = -mdot;

xdot(5) = xdot(5)*((xdot(5)>0) || ~withinPad);


xdot = xdot';

